//--------------------------------------------------------------------//
// Do-file: 04_Figures 
//
// Paper: "Little evidence that military policing reduces crime 
// or improves human security"
//
// Authors: Robert Blair and Michael Weintraub
//
// Date created: 2022-06-30
// Date last modified:2023-01-06
// 
// Notes:
// (1) This code builds the figures included in the main section and the appendix
// of the paper. Each figures includes the randomization inference p-value. To 
// calculate it, each regression is run 10,000 times with different simulations 
// of the randomization process. 
//
//--------------------------------------------------------------------//


					//////////////////////////////////////
					///								   ///
					///		Figures of Main Results	   ///
					///								   ///
					//////////////////////////////////////

					
//---------//
// -- Preliminaries 
//---------//

clear all
set more off
graph set window fontface "Times New Roman"


* Defining main computer path
global path "~/Dropbox/Asistencia de investigacion/Proyecto Cali/Mapas Cali/08_Data&Analysis/02_dofiles/v20220630_PaperReplication" 

		
* Set working directory
cd "${path}" 

	
* Uploading housekeeping
include "${path}/Code/01_Housekeeping.do"
	
		
* Increase max number of variables 
set maxvar 20000


			//---------------------------------------------------------//
			//-- Figure 2. Randomized military patrols had little to --//
			//		no effect on crime even while soldiers were		   //
			// 			physically present on the streets 	           //
			//---------------------------------------------------------//

//----//								 
//- Regressions with pre-strike data (during intervention)	
//----//
			
	use "$data/raw/admin_data_during.dta", clear
	
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male
	
	
//-- Column 1 (all crime)
	
	reg unw_crime2_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_unw_crime2_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store col1
	
	
//-- Column 2 (weekend crime)
	
	reg crime2_wend_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_crime2_wend_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store col2
	

//-- Column 3 (weekday crime)
	
	reg crime2_wday_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_crime2_wday_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store col3

	
//-- Column 4 (daytime crime)
	
	reg crime2_day_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_crime2_day_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store col4
		
		
//-- Column 5 (nighttime crime)
	
	reg crime2_night_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_crime2_night_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store col5
		
		
//-- Coefplot of all regressions

matrix T = J(5,3,.)
matrix S = J(5,3,.)

foreach n of numlist 1/5{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore col1
	} 
	
	if `n' == 2 {
		estimates restore col2
	} 
	
	if `n' == 3 {
		estimates restore col3
	}
	
	if `n' == 4 {
		estimates restore col4
	}
	
	if `n' == 5 {
		estimates restore col5
	}
	
	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "All crime" "Weekend crime" "Weekday crime" "Daytime crime" "Nighttime crime" 
	matrix colname T = coef ll ul 
	matrix rowname T = "All crime" "Weekend crime" "Weekday crime" "Daytime crime" "Nighttime crime" 

	
//-------//
//-- Run regressions for randomization inference (admin data pre-strike)
//-------//

* Define globals for regressions
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male 


* Run regressions

foreach var in unw_crime2_num crime2_wend_num crime2_wday_num crime2_day_num crime2_night_num {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			
	 
	//- Merge with outcomes, weights and controls
		
			* Merge with endline database
			merge 1:m manzana using "$data/raw/admin_data_during.dta"
				drop _merge
	
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri /// 
			i.barrio_code $geovars $blockdemovars cum_all_`var' [pweight=iweight]
		
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figure2_`var'.dta", replace
	
}
	

	
//------//					 
//-- Store p-values from randomization inference and export table 
//------//

loc n = 0	
foreach var in unw_crime2_num crime2_wend_num crime2_wday_num crime2_day_num crime2_night_num {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	use "$data/rand_inf/RI_figure2_`var'.dta", clear
	loc ++n
	
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore col1
	} 
	
	if `n' == 2 {
		estimates restore col2
	} 
	
	if `n' == 3 {
		estimates restore col3
	}
	
	if `n' == 4 {
		estimates restore col4
	}
	
	if `n' == 5 {
		estimates restore col5
	}
	

	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f")
	di `p_spill`n''
} 


//-- Export table 

	loc var1 `"\textbf{\shortstack{All \\ crime}}"'
	loc var2 `"\textbf{\shortstack{Weekend \\ crime}}"'
	loc var3 `"\textbf{\shortstack{Weekday \\ crime}}"'
	loc var4 `"\textbf{\shortstack{Daytime \\ crime}}"'
	loc var5 `"\textbf{\shortstack{Nighttime \\ crime}}"'
	
	esttab col1 col2 col3 col4 col5 ///
		using "$tables/table_figure_2.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${demovars}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${geovars}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var3'"' `"`var4'"' `"`var5'"')  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{5}{c}{\textbf{Admin data}} \\ 		///
			\cline{2-6} \addlinespace					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )	
	   
	
//------//					 
//-- Coefplot of all regressions
//------//	

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.07) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.07) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))	

graph export "$figures/figure_2.pdf", replace	
			
			
						 
			 //-------------------------------------------------------//
			 //-- 	Figure 3. Randomized military patrols did not   --//
			 //		improve and may have diminished perceptions of    //
			 // 		safety except among business owners 	      //
			 //-------------------------------------------------------//
						
//- Load data (endline)
					
	use "$data/raw/survey_endline.dta", clear

	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
		

//-- Column 1 (all security)

	reg i_securityallindex_std  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
		vce(cluster manzana_code) baselevels
	
	estimates store fig3_1
	estadd matrix table = r(table)
	
	qui tab treatment if e(sample), matcell(M)
	estadd scalar n_treat = M[2,1], replace 
	estadd scalar n_spill = M[3,1], replace 
	
	
//-- Column 2 (security residents)

	reg i_securityallindex_ps  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
		vce(cluster manzana_code) baselevels
		
	estimates store fig3_2
	estadd matrix table = r(table)
	
	qui tab treatment if e(sample), matcell(M)
	estadd scalar n_treat = M[2,1], replace 
	estadd scalar n_spill = M[3,1], replace 
		
		
//-- Column 3 (business owners)

	reg i_business  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
		vce(cluster manzana_code) baselevels
		
	estimates store fig3_3
	estadd matrix table = r(table)

	qui tab treatment if e(sample), matcell(M)
	estadd scalar n_treat = M[2,1], replace 
	estadd scalar n_spill = M[3,1], replace
	
		
//-- Coefplot of all regressions

	matrix T = J(4,3,.)
	matrix S = J(4,3,.)

	foreach n of numlist 1/3{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore fig3_1
	} 
	
	if `n' == 2 {
		estimates restore fig3_2
	} 
	
	if `n' == 3 {
		estimates restore fig3_3
	}
	
	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "All safety" "Personal safety" "Business safety" 
	matrix colname T = coef ll ul 
	matrix rowname T = "All safety" "Personal safety" "Business safety" 



//-------//
//-- Run regressions for randomization inference
//-------//

* Define globals for regressions
	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 


* Run regressions

foreach var in i_securityallindex_std i_securityallindex_ps i_businessindex_std {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
			
	 
	//- Merge with outcomes, weights and controls
		
			* Merge with endline database
			merge 1:m manzana_code using "$data/raw/survey_endline.dta"
				drop _merge
			
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figure3_`var'.dta", replace
	
}



//------//					 
//-- Store p-values from randomization inference and make table 
//------//

loc n = 0	
foreach var in i_securityallindex_std i_securityallindex_ps i_businessindex_std {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	use "$data/rand_inf/RI_figure3_`var'.dta", clear
	loc ++n
	
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore fig3_1
	} 
	
	if `n' == 2 {
		estimates restore fig3_2
	} 
	
	if `n' == 3 {
		estimates restore fig3_3
	}
	

	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f")
	di `p_spill`n''
} 


//-- Export table 

	loc var1 `"\textbf{\shortstack{All \\ safety}}"'
	loc var2 `"\textbf{\shortstack{Personal \\ safety}}"'
	loc var3 `"\textbf{\shortstack{Business \\ safety}}"'
	
	esttab fig3_1 fig3_2 fig3_3 ///
		using "$tables/table_figure_3.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${demovars}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${geovars}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover n_treat n_spill, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)" "Obs treatment" "Obs spillover")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var3'"')  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{3}{c}{\textbf{Survey data}} \\ 		///
			\cline{2-4} \addlinespace					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )		


	
//------//					 
//-- Coefplot of all regressions
//------//	

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.06) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.06) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
		
graph export "$figures/figure_3.pdf", replace	

		
	
				//---------------------------------------------------//
				//--   Figure 4. Randomized military patrols may   --//
				//  	 have exacerbated human rights abuses,       // 
				//			especially by the police 		 		 //
				//---------------------------------------------------//

//----//								 
//- Regressions with monitoring data	
//----//
				
	use "$data/raw/survey_monitoring.dta", clear
	
	global indiv_control age gender 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	
	
//-- Column 1 (abuses by the police)	
					
	reg abuse_police  i.treatment i.barrio_code ${geovars} ${indiv_control} [pweight=iweight], ///
		vce(cluster manzana_code) baselevels
	
	estimates store t3_1,  title("Police")
	
	estadd matrix table = r(table)
	
	
//-- Column 2 (abuses by the military)

	reg abuse_military  i.treatment i.barrio_code ${geovars} ${indiv_control} [pweight=iweight], ///
		vce(cluster manzana_code) baselevels
	
	estimates store t3_2,  title("Military")
	
	estadd matrix table = r(table)
	

//----//								 
//- Regressions with endline data	
//----//

	use "$data/raw/survey_endline.dta", clear

	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
		
	
//-- Column 3 (abuses by the police)

	reg abuse_police  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
		vce(cluster manzana_code) baselevels
		
	estimates store t3_3,  title("Police")
	
	estadd matrix table = r(table)
		
		
//-- Column 4 (military abuse)

	reg abuse_military  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels
		
	estimates store t3_4,  title("Military")
	
	estadd matrix table = r(table)
	
		
//-- Coefplot of all regressions

	matrix T = J(4,4,.)
	matrix S = J(4,4,.)

	foreach n of numlist 1/4{
		
		* Restore estimates of each regression
		
		if `n' == 1 {
			estimates restore t3_1
		} 
		
		if `n' == 2 {
			estimates restore t3_2
		} 
		
		if `n' == 3 {
			estimates restore t3_3
		}
		
		if `n' == 4 {
			estimates restore t3_4
		}
		
		* Store estimates in matrix 
		
		matrix t_table = e(table)
				matrix T[`n',1] = t_table[1,"1.treatment"]
				matrix T[`n',2] = t_table[5,"1.treatment"]
				matrix T[`n',3] = t_table[6,"1.treatment"]
				
				matrix S[`n',1] = t_table[1,"2.treatment"]
				matrix S[`n',2] = t_table[5,"2.treatment"]
				matrix S[`n',3] = t_table[6,"2.treatment"]
	}


	* Name matrix elements 
	
	matrix colname S = coef ll ul pval
	matrix rowname S = "Police abuse (monitoring)" "Military abuse (monitoring)" "Police abuse (survey)" "Military abuse (survey)"
	matrix colname T = coef ll ul pval
	matrix rowname T = "Police abuse (monitoring)" "Military abuse (monitoring)" "Police abuse (survey)" "Military abuse (survey)"
						
						
//-------//
//-- Run regressions for randomization inference (monitoring data)
//-------//

* Define globals for regressions
	global indiv_control age gender 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	

* Run regressions

foreach var in abuse_police abuse_military {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
	 
	//- Merge with outcomes
			merge 1:m manzana_code using "$data/raw/survey_monitoring.dta"
				keep if _merge == 3
				drop _merge
			
		
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${geovars} ${indiv_control} [pweight=iweight], ///
				vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figure4_`var'.dta", replace
	
}

						

//-------//
//-- Run regressions for randomization inference (endline data)
//-------//

* Define globals for regressions
	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 


* Run regressions

foreach var in abuse_police_end abuse_military_end {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
			
	 
	//- Merge with outcomes
			merge 1:m manzana_code using "$data/raw/survey_endline.dta"
				keep if _merge == 3
				drop _merge
		
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels
	
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figure4_`var'.dta", replace
	
}



//------//					 
//-- Store p-values from randomization inference and export the table 
//------//

loc n = 0	
foreach var in abuse_police abuse_military abuse_police_end abuse_military_end {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	use "$data/rand_inf/RI_figure4_`var'.dta", clear
	loc ++n
	
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore t3_1
	} 
	
	if `n' == 2 {
		estimates restore t3_2
	} 
	
	if `n' == 3 {
		estimates restore t3_3
	}
	
	if `n' == 4 {
		estimates restore t3_4
	}
	

	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f")
	di `p_spill`n''
} 


//-- Export table 

	loc var1 `"\textbf{\shortstack{Police \\ abuse}}"'
	loc var2 `"\textbf{\shortstack{Military \\ abuse}}"'
	loc var3 `"\textbf{\shortstack{Police \\ abuse}}"'
	loc var4 `"\textbf{\shortstack{Military \\ abuse}}"'
	
	esttab t3_1 t3_2 t3_3 t3_4 ///
		using "$tables/table_figure_4.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${demovars}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${geovars}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var3'"'  `"`var4'"')  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{2}{c}{\textbf{Monitoring survey}} & \multicolumn{2}{c}{\textbf{Endline survey}} \\  ///
			\cline{2-5} \addlinespace					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )		

	   
//------//
//-- Coefplot of all regressions
//------//

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.08) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.08) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 coeflabels("Police abuse (monitoring)" = "Police abuse" "Police abuse (survey)" = "Police abuse" /// 
						"Military abuse (monitoring)" = "Military abuse" "Military abuse (survey)" = "Military abuse") ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("Police abuse (monitoring)" "Military abuse (monitoring)" = `""{bf:During intervention}" "{bf:(monitoring survey)}""' ///
					  "Police abuse (survey)" "Military abuse (survey)" = `""{bf:After intervention}" "{bf:(endline survey)}""') ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))

graph export "$figures/figure_4.pdf", replace	
					  

	
					///////////////////////////////////////////
					///								   		///
					///		Figures for the appendix 		///
					///				of the paper			///
					///								   		///
					///////////////////////////////////////////

				//----------------------------------------------//
				//--  Figure A.6. Treatment effects on crime  --//
				// 	  	using weighted standardized index       //
				//----------------------------------------------//

//----//								 
//- Regressions with pre-strike data (during intervention)	
//----//
			
	use "$data/raw/admin_data_during.dta", clear
	
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male

	
//-- Regressions 
	
	reg w_crime2_num_std  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_w_crime2_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store ACP_noInt,  title("During intervention")
	

//----//								 
//- Regressions with post-strike data	
//----//	

	use "$data/raw/admin_data_after.dta", clear
	
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male

	
//-- Regressions 

	reg w_crime_num_std  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_w_crime_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store AdminShort2020,  title("After intervention")
	

//-----//		
//-- Coefplot of all regressions
//-----//
 
//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(2,3,.)
matrix S = J(2,3,.)

foreach n of numlist 1/2{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore ACP_noInt
	} 
	
	if `n' == 2 {
		estimates restore AdminShort2020
	} 

	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "during" "after"   
	matrix colname T = coef ll ul
	matrix rowname T = "during" "after"   

						
//-------//
//-- Run regressions for randomization inference (admin data pre-strike)
//-------//

* Define globals for regressions
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male 


* Run regressions

foreach var in w_crime2_num_std  {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			
	 
	//- Merge with outcomes
			merge 1:m manzana using "$data/raw/admin_data_during.dta"
				drop _merge
	
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri /// 
			i.barrio_code $geovars $blockdemovars cum_all_w_crime2_num [pweight=iweight]
		
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA6_`var'.dta", replace
	
}
			
			
//-------//
//-- Run regressions for randomization inference (admin data)
//-------//

* Define globals for regressions
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male


* Run regressions

foreach var in w_crime_num_std  {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			
	 
	//- Merge with outcomes
			merge 1:m manzana using "$data/raw/admin_data_after.dta"
				drop _merge
	
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri /// 
			i.barrio_code $geovars $blockdemovars cum_all_w_crime_num [pweight=iweight]
		
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA6_`var'.dta", replace
	
}
			


//------//					 
//-- Store p-values from randomization inference and export the table 
//------//

loc n = 0	
foreach var in w_crime2_num_std w_crime_num_std {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	use "$data/rand_inf/RI_figureA6_`var'.dta", clear
	loc ++n
	
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore ACP_noInt
	} 
	
	if `n' == 2 {
		estimates restore AdminShort2020
	} 
	

	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f")
	di `p_spill`n''
} 


//-- Export table 

	loc var1 `"\textbf{\shortstack{During \\intervention}}"'
	loc var2 `"\textbf{\shortstack{After \\intervention}}"'
	
	esttab ACP_noInt AdminShort2020  ///
		using "$tables/table_figure_A6.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${demovars}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${geovars}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"')  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{2}{c}{\textbf{Admin data}} \\  ///
			\cline{2-3} \addlinespace					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )	
	   
						
//-------//			
//-- Coefplot of all regressions
//-------//

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.06) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.06) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("during" "after" = `""{bf:Admin data}" "{bf:(Crime index)}""') ///
			 coeflabels("during" = "During intervention" "after" = "After intervention") ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
				
graph export "$figures/figure_A6.pdf", replace	
				
	
	
				//-----------------------------------------------------//
				//-- 	Figure A7. Treatment effects on crime by day --//
				//		and time using weighted standardized index     //
				//-----------------------------------------------------//

//----//								 
//- Regressions with pre-strike data (during intervention)	
//----//
			
	use "$data/raw/admin_data_during.dta", clear
	
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male
	
	
//-- Column 1 (all crime)
	
	reg w_crime2_num_std  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_w_crime2_num_std [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store col1
	
	
//-- Column 2 (weekend crime)
	
	reg w_crime2_wend_num_std  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_w_crime2_wend_num_std [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store col2
	

//-- Column 3 (weekday crime)
	
	reg w_crime2_wday_num_std  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_w_crime2_wday_num_std [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store col3

	
//-- Column 4 (daytime crime)
	
	reg w_crime2_day_num_std  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_w_crime2_day_num_std [pweight=iweight]
	
	
	estadd matrix table = r(table)
	estimates store col4
		
		
//-- Column 5 (nighttime crime)
	
	reg w_crime2_night_num_std  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_w_crime2_night_num_std [pweight=iweight]
		
	estadd matrix table = r(table)
	estimates store col5
		
		
//-- Coefplot of all regressions

matrix T = J(5,3,.)
matrix S = J(5,3,.)

foreach n of numlist 1/5{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore col1
	} 
	
	if `n' == 2 {
		estimates restore col2
	} 
	
	if `n' == 3 {
		estimates restore col3
	}
	
	if `n' == 4 {
		estimates restore col4
	}
	
	if `n' == 5 {
		estimates restore col5
	}
	
	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "All crime" "Weekend crime" "Weekday crime" "Daytime crime" "Nighttime crime" 
	matrix colname T = coef ll ul 
	matrix rowname T = "All crime" "Weekend crime" "Weekday crime" "Daytime crime" "Nighttime crime" 

	
//-------//
//-- Run regressions for randomization inference (during intervention)
//-------//

* Define globals for regressions
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male 


* Run regressions

foreach var in w_crime2_num_std w_crime2_wend_num_std w_crime2_wday_num_std /*
	*/ w_crime2_day_num_std w_crime2_night_num_std {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			
	 
	//- Merge with outcomes
			merge 1:m manzana using "$data/raw/admin_data_during.dta"
				drop _merge
	
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri /// 
			i.barrio_code $geovars $blockdemovars cum_all_`var' [pweight=iweight]
		
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA7_`var'.dta", replace
	
}


//------//					 
//-- Store p-values from randomization inference and export table 
//------//

loc n = 0	
foreach var in w_crime2_num_std w_crime2_wend_num_std w_crime2_wday_num_std w_crime2_day_num_std w_crime2_night_num_std  {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	use "$data/rand_inf/RI_figureA7_`var'.dta", clear
	loc ++n
	
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore col1
	} 
	
	if `n' == 2 {
		estimates restore col2
	} 
	
	if `n' == 3 {
		estimates restore col3
	}
	
	if `n' == 4 {
		estimates restore col4
	}
	
	if `n' == 5 {
		estimates restore col5
	}
	

	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f")
	di `p_spill`n''
} 


//-- Export table 

	loc var1 `"\textbf{\shortstack{All \\ crime}}"'
	loc var2 `"\textbf{\shortstack{Weekend \\ crime}}"'
	loc var3 `"\textbf{\shortstack{Weekday \\ crime}}"'
	loc var4 `"\textbf{\shortstack{Daytime \\ crime}}"'
	loc var5 `"\textbf{\shortstack{Nighttime \\ crime}}"'
	
	esttab col1 col2 col3 col4 col5 ///
		using "$tables/table_figure_A7.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${demovars}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${geovars}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var3'"' `"`var4'"' `"`var5'"')  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{5}{c}{\textbf{Admin data}} \\ 		///
			\cline{2-6} \addlinespace					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )	
	   
	
//------//					 
//-- Coefplot of all regressions
//------//	

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.07) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.07) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) xtitle("ITT effect size", size(small)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) ///
			 ylabel(, nogrid) ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))	

graph export "$figures/figure_A7.pdf", replace	
			
			
	
				//--------------------------------------------------//
				//--  Figure A8. Treatment effects on violent and   //
				// 		  non-violent crime in admin data   		//
				//--------------------------------------------------//

//----//								 
//- Regressions with pre-strike data (during intervention)	
//----//
			
	use "$data/raw/admin_data_during.dta", clear
	
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male

	
//-- Regression column 1
	
	reg viol2_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_viol2_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store viol2_dur,  title("Violent crime")
	
	
//-- Regression column 2
	
	reg nonviol2_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_nonviol2_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store nonviol2_dur,  title("Non-Violent crime")
		
	
//----//								 
//- Regressions with post-strike data (after intervention)
//----//	

	use "$data/raw/admin_data_after.dta", clear 
	
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male

	
//-- Regression column 3
	
	reg viol_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_viol_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store AdminShort2020_viol,  title("Violent crime")
	
	
//-- Regression column 4

	reg nonviol_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_nonviol_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store AdminShort2020_nonviol,  title("Non-Violent crime")
	
							
//-----//		
//-- Coefplot of all regressions
//-----//
 
//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(4,3,.)
matrix S = J(4,3,.)

foreach n of numlist 1/4{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore viol2_dur
	} 
	
	if `n' == 2 {
		estimates restore nonviol2_dur
	} 
	
	if `n' == 3 {
		estimates restore AdminShort2020_viol
	} 
	
	if `n' == 4 {
		estimates restore AdminShort2020_nonviol
	} 

	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "viol_dur" "noviol_dur" "viol_after" "noviol_after"   
	matrix colname T = coef ll ul
	matrix rowname T = "viol_dur" "noviol_dur" "viol_after" "noviol_after" 
						
						
//-------//
//-- Run regressions for randomization inference (admin data pre-strike)
//-------//

* Define globals for regressions
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male 


* Run regressions

foreach var in viol2_num nonviol2_num {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			
	 
	//- Merge with outcomes, weights and controls
		
			* Merge with endline database
			merge 1:m manzana using "$data/raw/admin_data_during.dta"
				drop _merge
	
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri /// 
			i.barrio_code $geovars $blockdemovars cum_all_`var' [pweight=iweight]
		
		
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA8_`var'.dta", replace
	
}
						
						
						
//-------//
//-- Run regressions for randomization inference (admin data)
//-------//

* Define globals for regressions
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male


* Run regressions

foreach var in viol_num nonviol_num  {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			
	 
	//- Merge with outcomes
			merge 1:m manzana using "$data/raw/admin_data_after.dta"
				drop _merge
	
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri /// 
			i.barrio_code $geovars $blockdemovars cum_all_`var' [pweight=iweight]
		
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA8_`var'.dta", replace
	
}
							

//------//					 
//-- Store p-values from randomization inference and export the table 
//------//

loc n = 0	
foreach var in viol2_num nonviol2_num viol_num nonviol_num  {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	use "$data/rand_inf/RI_figureA8_`var'.dta", clear
	loc ++n
	
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore viol2_dur
	} 
	
	if `n' == 2 {
		estimates restore nonviol2_dur
	} 
	
	if `n' == 3 {
		estimates restore AdminShort2020_viol
	} 
	
	if `n' == 4 {
		estimates restore AdminShort2020_nonviol
	} 
	

	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f")
	di `p_spill`n''
} 


//-- Export table 

	loc var1 `"\textbf{\shortstack{Violent \\ crime}}"'
	loc var2 `"\textbf{\shortstack{Non-violent \\ crime}}"'
	loc var3 `"\textbf{\shortstack{Violent \\ crime}}"'
	loc var4 `"\textbf{\shortstack{Non-violent \\ crime}}"'
	
	esttab viol2_dur nonviol2_dur AdminShort2020_viol AdminShort2020_nonviol ///
		using "$tables/table_figure_A8.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${demovars}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${geovars}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var3'"' `"`var4'"')  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{2}{c}{\textbf{During intervention}} & \multicolumn{2}{c}{\textbf{After intervention}} \\ 		///
			\cline{2-5} \addlinespace					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )	
	   							

//------//					 
//-- Coefplot of all regressions
//------//

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.08) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.08) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) xtitle("ITT effect size", size(small)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) ///
			 ylabel(, nogrid) ///
			 groups("viol_dur" "noviol_dur" = `""{bf:During intervention}" "{bf:(admin data)}""' ///
					"viol_after" "noviol_after" = `""{bf:After intervention}" "{bf:(admin data)}""') ///
			 coeflabels("viol_dur" = "Violent crime" "noviol_dur" = "Non-violent crime" /// 
						"viol_after" = "Violent crime" "noviol_after" = "Non-violent crime") ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
								
graph export "$figures/figure_A8.pdf", replace	
									
						
						
					//------------------------------------------------//
					//--  Figure A9. Treatment effects on crime in -- //
					// 	   administrative data disaggregated by type  //
					//------------------------------------------------//
		
//----//								 
//- Regressions with pre-strike data (during intervention)	
//----//
			
	use "$data/raw/admin_data_during.dta", clear
	
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male

	
//-- Regression column 1
	
	reg armed_robberies_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_robberies_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store during2019_armrobb
	
	
//-- Regression column 2
	
	reg thefts_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_robberies_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store during2019_thefts
	
	
//-- Regression column 3
	
	reg homicides_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_homicides_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store during2019_homi
	
	
//-- Regression column 4
	
	reg drugdeal_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_drugdeal_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store during2019_drug
	
	
//-- Regression column 5
	
	reg wbearing_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_wbearing_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store during2019_weap
				
		
	
//----//								 
//- Regressions with post-strike data (after intervention)
//----//	

	use "$data/raw/admin_data_after.dta", clear
	
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male

	
//-- Regression column 6

	reg armed_robberies_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_robberies_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store AdminShort2020_armrobb,  title("Violent crime")
	
	
//-- Regression column 7

	reg thefts_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_robberies_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store AdminShort2020_thefts,  title("Violent crime")
	
	
//-- Regression column 8

	reg homicides_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_homicides_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store AdminShort2020_homi,  title("Non-Violent crime")
	
	
//-- Regression column 9

	reg drugdeal_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_drugdeal_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store AdminShort2020_drug,  title("Violent crime")
	
	
//-- Regression column 10

	reg wbearing_num  i.treatment /// 
		i.barrio_code $geovars $blockdemovars cum_all_wbearing_num [pweight=iweight]
	
	estadd matrix table = r(table)
	estimates store AdminShort2020_weap,  title("Non-Violent crime")
		

//-----//		
//-- Coefplot of all regressions
//-----//
 
//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(10,3,.)
matrix S = J(10,3,.)

foreach n of numlist 1/10{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore during2019_armrobb
	} 
	
	if `n' == 2 {
		estimates restore during2019_thefts
	} 
	
	if `n' == 3 {
		estimates restore during2019_homi
	} 
	
	if `n' == 4 {
		estimates restore during2019_drug
	} 
	
	if `n' == 5 {
		estimates restore during2019_weap
	} 
	
	if `n' == 6 {
		estimates restore AdminShort2020_armrobb
	} 
	
	if `n' == 7 {
		estimates restore AdminShort2020_thefts
	} 
	if `n' == 8 {
		estimates restore AdminShort2020_homi
	} 
	
	if `n' == 9 {
		estimates restore AdminShort2020_drug
	} 
	
	if `n' == 10 {
		estimates restore AdminShort2020_weap
	} 


	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "arm_rob_dur" "theft_dur" "homi_dur" "drug_dur" "weap_dur"  "arm_rob_aft" "theft_aft" "homi_aft" "drug_aft" "weap_aft"  
	matrix colname T = coef ll ul
	matrix rowname T = "arm_rob_dur" "theft_dur" "homi_dur" "drug_dur" "weap_dur"  "arm_rob_aft" "theft_aft" "homi_aft" "drug_aft" "weap_aft"

	
		
//-------//
//-- Run regressions for randomization inference (admin data pre-strike)
//-------//

* Define globals for regressions
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male 
	
	
* Run regressions

foreach var in armed_robberies_num thefts_num homicides_num drugdeal_num wbearing_num {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			
	 
	//- Merge with outcomes
			merge 1:m manzana using "$data/raw/admin_data_during.dta"
				drop _merge
	
	
	//-- Define cummulative crime var
	
		if ("`var'" == "armed_robberies_num" | "`var'" == "thefts_num") {
			local var_cum robberies_num
			} 
		
		else {
			local var_cum `var'
			}
	
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri /// 
			i.barrio_code $geovars $blockdemovars cum_all_`var_cum' [pweight=iweight]
		
		
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA9_`var'.dta", replace
	
}



//-------//
//-- Run regressions for randomization inference (admin data poststrike)
//-------//

* Define globals for regressions
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	global blockdemovars block_age block_educ block_pct_male


* Run regressions

foreach var in armed_robberies_num thefts_num homicides_num drugdeal_num wbearing_num {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			
	 
	//- Merge with outcomes
			merge 1:m manzana using "$data/raw/admin_data_after.dta"
				drop _merge
			
			
	//-- Define cummulative crime var
	
		if ("`var'" == "armed_robberies_num" | "`var'" == "thefts_num") {
			local var_cum robberies_num
			} 
		
		else {
			local var_cum `var'
			}
			
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri /// 
			i.barrio_code $geovars $blockdemovars cum_all_`var_cum' [pweight=iweight]
		
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA9_after_`var'.dta", replace
	
}
	

	
//------//					 
//-- Store p-values from randomization inference and export the table 
//------//

loc n = 0	
foreach var in armed_robberies_num thefts_num homicides_num drugdeal_num wbearing_num ///
			   armed_robberies_num thefts_num homicides_num drugdeal_num wbearing_num {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	loc ++n
	
	
		if (`n' <= 5) {
			use "$data/rand_inf/RI_figureA9_`var'.dta", clear
			} 
		
		else {
			use "$data/rand_inf/RI_figureA9_after_`var'.dta", clear
			}
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore during2019_armrobb
	} 
	
	if `n' == 2 {
		estimates restore during2019_thefts
	} 
	
	if `n' == 3 {
		estimates restore during2019_homi
	} 
	
	if `n' == 4 {
		estimates restore during2019_drug
	} 
	
	if `n' == 5 {
		estimates restore during2019_weap
	} 
	
	if `n' == 6 {
		estimates restore AdminShort2020_armrobb
	} 
	
	if `n' == 7 {
		estimates restore AdminShort2020_thefts
	} 
	if `n' == 8 {
		estimates restore AdminShort2020_homi
	} 
	
	if `n' == 9 {
		estimates restore AdminShort2020_drug
	} 
	
	if `n' == 10 {
		estimates restore AdminShort2020_weap
	} 
	

	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f")
	di `p_spill`n''
} 
	

//-- Export table 

	loc var1 `"\textbf{\shortstack{Armed \\robbery}}"'
	loc var2 `"\textbf{\shortstack{Thefts}}"'
	loc var3 `"\textbf{\shortstack{Homicide}}"'
	loc var4 `"\textbf{\shortstack{Drug possession}}"'
	loc var5 `"\textbf{\shortstack{Weapons possession}}"'
	loc var6 `"\textbf{\shortstack{Armed \\robbery}}"'
	loc var7 `"\textbf{\shortstack{Thefts}}"'
	loc var8 `"\textbf{\shortstack{Homicide}}"'
	loc var9 `"\textbf{\shortstack{Drug possession}}"'
	loc var10 `"\textbf{\shortstack{Weapons possession}}"'
	
	esttab during2019_armrobb during2019_thefts during2019_homi during2019_drug during2019_weap ///
		AdminShort2020_armrobb AdminShort2020_thefts AdminShort2020_homi AdminShort2020_drug AdminShort2020_weap ///
		using "$tables/table_figure_A9.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${demovars}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${geovars}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var3'"' `"`var4'"' `"`var5'"' `"`var6'"' `"`var7'"' `"`var8'"' `"`var9'"' `"`var10'"' )  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{5}{c}{\textbf{During intervention}} & \multicolumn{5}{c}{\textbf{After intervention}} \\ 		///
			\cline{2-6} \cline{7-11} \addlinespace					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )	
	   		
	
//-------//
//-- Coefplot of all regressions
//-------//

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.136) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.136) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("arm_rob_dur" "weap_dur" = `""{bf:During intervention}" "{bf:(admin data)}""' ///
					"arm_rob_aft" "weap_aft" = `""{bf:After intervention}" "{bf:(admin data)}""') ///
			 coeflabels("arm_rob_dur" = "Armed robbery" "theft_dur" = "Theft" "homi_dur" = "Homicide" /// 
						"drug_dur" = `""Drug" "dealing""' "weap_dur" = `""Possession of" "firearms""' /// 
						"arm_rob_aft" = "Armed robbery" "theft_aft" = "Theft" "homi_aft" = "Homicide" /// 
						"drug_aft" = `""Drug" "dealing""'  "weap_aft" = `""Possession of" "firearms""') ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
		
graph export "$figures/figure_A9.pdf", replace	
			
		
		
			//------------------------------------------------//
			//--	 	   	Figure A10. Part1 		    	--//
			//  Treatment effects on crime victimization in   //
			//  	survey data disaggregated by type         //
			//------------------------------------------------//

//-------//
//-- Run regressions for randomization inference (endline data)
//-------//

// Note: the code below runs the randomization inference for all variables 
// of parts 1 and 2 of figure A10. Further down separate figures are produced. 

* Define globals for regressions
	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 


* Run regressions

foreach var in vandalism_vicd_during2 robberies_vicd_during2 homicides_vicd_during2 ///
			homicides_intents_vicd_during2 assault_vicd_during2 car_robberie_vicd_during2 ///
			house_robberie_vicd_during2 gangs_vicd_during2 extorsion_vicd_during2 ///
			vandalism_vicd_after2 robberies_vicd_after2 homicides_vicd_after2 /// 
			homicides_intents_vicd_after2 assault_vicd_after2 car_robberie_vicd_after2 ///
			house_robberie_vicd_after2 gangs_vicd_after2 extorsion_vicd_after2 {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
			
	 
	//- Merge with outcomes
			merge 1:m manzana_code using "$data/raw/survey_endline.dta"
				drop _merge
	
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA10_`var'.dta", replace
	
}


							//------------------------//
							//-- Figure A10. Part1 	--//
							//------------------------//
						
//----//								 
//- Regressions with endline data	
//----//
						
	use "$data/raw/survey_endline.dta", clear

	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 

	
//-- Column 1
	
	reg vandalism_vicd_during2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_vandalism
			

//-- Column 2
	
	reg robberies_vicd_during2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_robberies

	
//-- Column 3
	
	reg homicides_vicd_during2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_homi

	
//-- Column 4
	
	reg homicides_intents_vicd_during2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_homint

	
//-- Column 5
	
	reg assault_vicd_during2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_assault
	
					 
//-- Column 6
	
	reg car_robberie_vicd_during2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_crobberie


//-- Column 7
	
	reg house_robberie_vicd_during2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_hrobberie
	
	
//-- Column 8
	
	reg gangs_vicd_during2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_gangs
					 
					 
//-- Column 9
	
	reg extorsion_vicd_during2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_extorsion
					 
						 
//-----//		
//-- Coefplot of all regressions
//-----//
 
//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(9,3,.)
matrix S = J(9,3,.)

foreach n of numlist 1/9{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore c_dur2_vandalism
	} 
	
	if `n' == 2 {
		estimates restore c_dur2_robberies
	} 
	
	if `n' == 3 {
		estimates restore c_dur2_homi
	} 
	
	if `n' == 4 {
		estimates restore c_dur2_homint
	} 
	
	if `n' == 5 {
		estimates restore c_dur2_assault
	} 
	
	if `n' == 6 {
		estimates restore c_dur2_crobberie
	} 
	
	if `n' == 7 {
		estimates restore c_dur2_hrobberie
	} 
	
	if `n' == 8 {
		estimates restore c_dur2_gangs
	}
	
	if `n' == 9 {
		estimates restore c_dur2_extorsion
	}


	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "vand" "theft" "homi" "att_homi" "arm_rob" "mot_rob" "burg" "gang" "ext"
	matrix colname T = coef ll ul
	matrix rowname T = "vand" "theft" "homi" "att_homi" "arm_rob" "mot_rob" "burg" "gang" "ext"
						
						
//------//					 
//-- Store p-values from randomization inference (figure part 1)
//------//

loc n = 0	
foreach var in vandalism_vicd_during2 robberies_vicd_during2 homicides_vicd_during2 ///
			homicides_intents_vicd_during2 assault_vicd_during2 car_robberie_vicd_during2 ///
			house_robberie_vicd_during2 gangs_vicd_during2 extorsion_vicd_during2 {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	use "$data/rand_inf/RI_figureA10_`var'.dta", clear
	loc ++n
	
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore c_dur2_vandalism
	} 
	
	if `n' == 2 {
		estimates restore c_dur2_robberies
	} 
	
	if `n' == 3 {
		estimates restore c_dur2_homi
	} 
	
	if `n' == 4 {
		estimates restore c_dur2_homint
	} 
	
	if `n' == 5 {
		estimates restore c_dur2_assault
	} 
	
	if `n' == 6 {
		estimates restore c_dur2_crobberie
	} 
	
	if `n' == 7 {
		estimates restore c_dur2_hrobberie
	} 
	
	if `n' == 8 {
		estimates restore c_dur2_gangs
	}
	
	if `n' == 9 {
		estimates restore c_dur2_extorsion
	}
	

	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f")
	di `p_spill`n''
} 
				
				
//-- Export table (figure part 1)

	loc var1 `"\textbf{\shortstack{Vandalism}}"'
	loc var2 `"\textbf{\shortstack{Thefts}}"'
	loc var3 `"\textbf{\shortstack{Homicide}}"'
	loc var4 `"\textbf{\shortstack{Attempted \\homicides}}"'
	loc var5 `"\textbf{\shortstack{Armed \\robbery}}"'
	loc var6 `"\textbf{\shortstack{Motor vehicle theft}}"'
	loc var7 `"\textbf{\shortstack{Burlgary}}"'
	loc var8 `"\textbf{\shortstack{Gang activity}}"'
	loc var9 `"\textbf{\shortstack{Extortion}}"'
	
	esttab c_dur2_vandalism c_dur2_robberies c_dur2_homi c_dur2_homint c_dur2_assault ///
		c_dur2_crobberie c_dur2_hrobberie c_dur2_gangs c_dur2_extorsion ///
		using "$tables/table_figure_A10_part1.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${demovars}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${geovars}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var3'"' `"`var4'"' `"`var5'"' `"`var6'"' `"`var7'"' `"`var8'"' `"`var9'"' `"`var10'"' )  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{9}{c}{\textbf{During intervention}} \\ 		///
			\cline{2-10} \addlinespace					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )	
	   	
						
						
//-- Coefplot of all regressions (figure part 1)

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.13) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.13) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("vand" "ext" = `""{bf:Crime victimization }" "{bf: during intervention}" "{bf: (endline survey)}""' ) ///
			 coeflabels("vand" = "Vandalism" "theft" = "Theft" /// 
						"homi" = "Homicides" "att_homi" = `""Attempted" "homicides""' /// 
						"arm_rob" = `""Armed" "robbery""' "mot_rob" = `""Motor vehicle" "theft""' /// 
						"burg" = "Burglary"  "gang" = `""Gang" "activity""' ///
						"ext" = "Extortion") ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
	
graph export "$figures/figure_A10_part1.pdf", replace	
						

						
			//------------------------------------------------//
			//--	 		Figure A10. Part 2 				--//
			//  Treatment effects on crime victimization in   //
			//   survey data disaggregated by type (cont.)    //
			//------------------------------------------------//

//----//								 
//- Regressions with endline data
//----//
						
	use "$data/raw/survey_endline.dta", clear

	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 

	
//-- Column 1
	
	reg vandalism_vicd_after2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_vandalism
			

//-- Column 2
	
	reg robberies_vicd_after2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_robberies

	
//-- Column 3
	
	reg homicides_vicd_after2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_homi

	
//-- Column 4
	
	reg homicides_intents_vicd_after2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_homint

	
//-- Column 5
	
	reg assault_vicd_after2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_assault
	
					 
//-- Column 6
	
	reg car_robberie_vicd_after2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_crobberie


//-- Column 7
	
	reg house_robberie_vicd_after2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_hrobberie
	
	
//-- Column 8
	
	reg gangs_vicd_after2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_gangs
					 
					 
//-- Column 9
	
	reg extorsion_vicd_after2  i.treatment i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels

	estadd matrix table = r(table)
	estimates store c_dur2_extorsion
					 

						 
//-----//		
//-- Coefplot of all regressions
//-----//
 
//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(9,3,.)
matrix S = J(9,3,.)

foreach n of numlist 1/9{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore c_dur2_vandalism
	} 
	
	if `n' == 2 {
		estimates restore c_dur2_robberies
	} 
	
	if `n' == 3 {
		estimates restore c_dur2_homi
	} 
	
	if `n' == 4 {
		estimates restore c_dur2_homint
	} 
	
	if `n' == 5 {
		estimates restore c_dur2_assault
	} 
	
	if `n' == 6 {
		estimates restore c_dur2_crobberie
	} 
	
	if `n' == 7 {
		estimates restore c_dur2_hrobberie
	} 
	
	if `n' == 8 {
		estimates restore c_dur2_gangs
	}
	
	if `n' == 9 {
		estimates restore c_dur2_extorsion
	}


	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "vand" "theft" "homi" "att_homi" "arm_rob" "mot_rob" "burg" "gang" "ext"
	matrix colname T = coef ll ul
	matrix rowname T = "vand" "theft" "homi" "att_homi" "arm_rob" "mot_rob" "burg" "gang" "ext"

						
//------//					 
//-- Store p-values from randomization inference (figure part 2)
//------//

loc n = 0	
foreach var in vandalism_vicd_after2 robberies_vicd_after2 homicides_vicd_after2 /// 
				homicides_intents_vicd_after2 assault_vicd_after2 car_robberie_vicd_after2 ///
				house_robberie_vicd_after2 gangs_vicd_after2 extorsion_vicd_after2 {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	use "$data/rand_inf/RI_figureA10_`var'.dta", clear
	loc ++n
	
	
	* Restore estimates of each regression
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore c_dur2_vandalism
	} 
	
	if `n' == 2 {
		estimates restore c_dur2_robberies
	} 
	
	if `n' == 3 {
		estimates restore c_dur2_homi
	} 
	
	if `n' == 4 {
		estimates restore c_dur2_homint
	} 
	
	if `n' == 5 {
		estimates restore c_dur2_assault
	} 
	
	if `n' == 6 {
		estimates restore c_dur2_crobberie
	} 
	
	if `n' == 7 {
		estimates restore c_dur2_hrobberie
	} 
	
	if `n' == 8 {
		estimates restore c_dur2_gangs
	}
	
	if `n' == 9 {
		estimates restore c_dur2_extorsion
	}

	

	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f")
	di `p_spill`n''
} 


//-- Export table (part 2)

	loc var1 `"\textbf{\shortstack{Vandalism}}"'
	loc var2 `"\textbf{\shortstack{Thefts}}"'
	loc var3 `"\textbf{\shortstack{Homicide}}"'
	loc var4 `"\textbf{\shortstack{Attempted \\homicides}}"'
	loc var5 `"\textbf{\shortstack{Armed \\robbery}}"'
	loc var6 `"\textbf{\shortstack{Motor vehicle theft}}"'
	loc var7 `"\textbf{\shortstack{Burglary}}"'
	loc var8 `"\textbf{\shortstack{Gang activity}}"'
	loc var9 `"\textbf{\shortstack{Extortion}}"'
	
	esttab c_dur2_vandalism c_dur2_robberies c_dur2_homi c_dur2_homint c_dur2_assault ///
		c_dur2_crobberie c_dur2_hrobberie c_dur2_gangs c_dur2_extorsion ///
		using "$tables/table_figure_A10_part2.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${demovars}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${geovars}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var3'"' `"`var4'"' `"`var5'"' `"`var6'"' `"`var7'"' `"`var8'"' `"`var9'"' `"`var10'"' )  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{9}{c}{\textbf{After intervention}} \\ 		///
			\cline{2-10} \addlinespace					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )	
	   							

//-- Coefplot of all regressions (figure part 2)

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.13) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.13) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("vand" "ext" = `""{bf:Crime victimization }" "{bf: after intervention}" "{bf: (endline survey)}""' ) ///
			 coeflabels("vand" = "Vandalism" "theft" = "Theft" /// 
						"homi" = "Homicides" "att_homi" = `""Attempted" "homicides""' /// 
						"arm_rob" = `""Armed" "robbery""' "mot_rob" = `""Motor vehicle" "theft""' /// 
						"burg" = "Burglary"  "gang" = `""Gang" "activity""' ///
						"ext" = "Extortion") ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
	
graph export "$figures/figure_A10_part2.pdf", replace	
						
												
		
			//------------------------------------------------//
			//--		Figure A11. Part 1 and 2 			--//
			//  Treatment effects on crime witnessing in      //
			//  	survey data disaggregated by type         //
			//------------------------------------------------//

// Note: the code below runs the randomization inference for all variables 
// of parts 1 and 2 of figure A10. Further down separate figures are produced.  

//-------//
//-- Run regressions for randomization inference (endline data)
//-------//

* Define globals for regressions
	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 


* Run regressions

foreach var in wit2_prostitu_std wit2_drugtrade_std wit2_drugconsump_std /*
							  */ wit2_alcohol_std wit2_vandalism_std wit2_robberies_std /*
							  */ wit2_homicide_std wit2_homicide_intents_std wit2_assault_std /*
							  */ wit2_car_robberie_std wit2_house_robberie_std wit2_domestic_vio_std /*
							  */ wit2_gangs_std wit2_iweaponcarrying_std wit2_extorsion_std {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
			
	 
	//- Merge with outcomes, weights and controls
		
			* Merge with endline database
			merge 1:m manzana_code using "$data/raw/survey_endline.dta"
				drop _merge
	
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA11_`var'.dta", replace
	
}
	
	
						//------------------------//
						//-- Figure A11. Part1 	--//
						//------------------------//
		
//----//								 
//- Regressions with endline data	
//----//
						
	use "$data/raw/survey_endline.dta", clear

	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	
	
//-- Run all regressions 

	global allcrime wit2_prostitu_std wit2_drugtrade_std wit2_drugconsump_std /*
							  */ wit2_alcohol_std wit2_vandalism_std wit2_robberies_std /*
							  */ wit2_homicide_std wit2_homicide_intents_std wit2_assault_std /*
							  */ wit2_car_robberie_std wit2_house_robberie_std wit2_domestic_vio_std /*
							  */ wit2_gangs_std wit2_iweaponcarrying_std wit2_extorsion_std
				
	loc j = 1
	foreach x of varlist $allcrime {
					
					di "Loop of variable `x' number `j'"
					
					 reg `x' i.treatment ///
								i.barrio_code  ${demovars} ${geovars}  [pweight=iweight], /// 
								vce(cluster manzana_code) baselevels
						
						estimates store col`j'
						estadd matrix table = r(table)
						
							loc ++j		
	}
	

//-----//		
//-- Coefplot of regressions 1-8
//-----//
 
//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(8,3,.)
matrix S = J(8,3,.)

foreach n of numlist 1/8{
	
	* Restore estimates of each regression
	
	estimates restore col`n'
	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "col1" "col2" "col3" "col4" "col5" "col6" "col7" "col8"
	matrix colname T = coef ll ul
	matrix rowname T = "col1" "col2" "col3" "col4" "col5" "col6" "col7" "col8"

							
//------//					 
//-- Store p-values from randomization inference (figure part 1)
//------//

loc n = 0	
foreach var in wit2_prostitu_std wit2_drugtrade_std wit2_drugconsump_std /*
							  */ wit2_alcohol_std wit2_vandalism_std wit2_robberies_std /*
							  */ wit2_homicide_std wit2_homicide_intents_std {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	use "$data/rand_inf/RI_figureA11_`var'.dta", clear
	loc ++n
	
	
	* Restore estimates of each regression
	
	estimates restore col`n'

	
	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f")
	di `p_spill`n''
} 
						

//-- Export table (figure part 1)
	
	loc var1 `"\textbf{\shortstack{Prostitution}}"'
	loc var2 `"\textbf{\shortstack{Drug dealing}}"'
	loc var3 `"\textbf{\shortstack{Drug consumption}}"'
	loc var4 `"\textbf{\shortstack{Alcohol consumption}}"'
	loc var5 `"\textbf{\shortstack{Vandalism}}"'
	loc var6 `"\textbf{\shortstack{Theft}}"'
	loc var7 `"\textbf{\shortstack{Homicide}}"'
	loc var8 `"\textbf{\shortstack{Attempted homicide}}"'
	
	esttab col1 col2 col3 col4 col5 col6 col7 col8 ///
		using "$tables/table_figure_A11_part1.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${demovars}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${geovars}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var3'"' `"`var4'"' `"`var5'"' `"`var6'"' `"`var7'"' `"`var8'"' )  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{8}{c}{\textbf{After intervention}} \\ 		///
			\cline{2-9} \addlinespace					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )
	   
	   
//-- Coefplot of all regressions (figure part 1)

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.13) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.13) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("col1" "col8" = `""{bf:Crime witnessing }" "{bf:after intervention}" "{bf: (endline survey)}""' ) ///
			 coeflabels("col1" = "Prostitution" "col2" = `""Drug" "dealing""' /// 
						"col3" = `""Drug" "consumption""' "col4" = `""Alcohol" "consumption""' ///
						"col5" = "Vandalism" "col6" = "Theft" "col7" = "Homicide" ///
						"col8" = `""Attempted" "homicide""')  ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
	
graph export "$figures/figure_A11_part1.pdf", replace	
			
			
			
					//------------------------//
					//-- Figure A11. Part2 	--//
					//------------------------//					
					
//-----//		
//-- Coefplot of regressions 9-15
//-----//
 
//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(7,3,.)
matrix S = J(7,3,.)

foreach j of numlist 9/15{
	
	* Restore estimates of each regression
	
	local n = `j'-8
	estimates restore col`j'
	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "col1" "col2" "col3" "col4" "col5" "col6" "col7" 
	matrix colname T = coef ll ul
	matrix rowname T = "col1" "col2" "col3" "col4" "col5" "col6" "col7" 
					
						
//------//					 
//-- Store p-values from randomization inference (figure part 2)
//------//
							  
loc n = 8	
foreach var in wit2_assault_std wit2_car_robberie_std wit2_house_robberie_std ///
			wit2_domestic_vio_std wit2_gangs_std wit2_iweaponcarrying_std /// 
			wit2_extorsion_std {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	use "$data/rand_inf/RI_figureA11_`var'.dta", clear
	loc ++n
	
	
	* Restore estimates of each regression
	
	estimates restore col`n'

	
	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f") 
	di `p_spill`n''
} 
						

	
//-- Export table (figure part 2)
	
	loc var1 `"\textbf{\shortstack{Armed \\robbery}}"'
	loc var2 `"\textbf{\shortstack{Motor vehicle theft}}"'
	loc var3 `"\textbf{\shortstack{Burglary}}"'
	loc var4 `"\textbf{\shortstack{Domestic violence}}"'
	loc var5 `"\textbf{\shortstack{Gang activity}}"'
	loc var6 `"\textbf{\shortstack{Possession of firearms}}"'
	loc var7 `"\textbf{\shortstack{Extortion}}"'
	
	esttab col9 col10 col11 col12 col13 col14 col15 ///
		using "$tables/table_figure_A11_part2.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${demovars}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${geovars}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var3'"' `"`var4'"' `"`var5'"' `"`var6'"' `"`var7'"' )  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{7}{c}{\textbf{After intervention}} \\ 		///
			\cline{2-8} \addlinespace					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )	
	   

//-- Coefplot of all regressions (figure part 2)

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.12) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.12) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("col1" "col7" = `""{bf:Crime witnessing }" "{bf:after intervention}" "{bf: (endline survey)}""' ) ///
			 coeflabels("col1" =  `""Armed" "robbery""' "col2" = `""Motor vehicle" "theft""' /// 
						"col3" = "Burglary" "col4" = `""Domestic" "violence""' ///
						"col5" = `""Gang" "activity""' "col6" = `""Possession" "of firearms""' ///
						"col7" = "Extortion") ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
	
graph export "$figures/figure_A11_part2.pdf", replace	
			
							
					
				//-------------------------------------------------//
				//--	Figure A23. Treatment effects on abuses  --//
				//			disaggregated by date of survey	       //
				//-------------------------------------------------//

//-- Load data on survey data 
	
use "$data/raw/survey_monitoring", clear
global indiv_control age gender 
global geovars number_buildings_sampling area bat_min cai_min ptr_min 	
		
//-----//
//-- Regressions with abuses during the intervention (before November 19)
//-----//

//-- Run regressions
 
		preserve
			keep if date < mdy(11,19,2019)
			
				global allabuses physical_abuse_police verbal_abuse_police 		 /*
							  */ physical_abuse_military verbal_abuse_military 
				
				loc j = 1
				foreach x of varlist $allabuses {
					
					if (regexm("`x'", "police") == 1) loc ti "Police"
					if (regexm("`x'", "milita") == 1) loc ti "Military"
			
					 reg `x' i.treatment ///
								i.barrio_code $geovars $indiv_control  [pweight=iweight], vce(cluster manzana_code)
								
						estimates store abud`j',  title("`ti'")
						estadd matrix table = r(table)
						
							loc ++j		
				}
			
		restore

			
			
//-----//
//-- Regressions with abuses after the intervention (after November 19)
//-----//

//-- Run regressions

		preserve
			keep if date >= mdy(11,19,2019)		
			
				global allabuses physical_abuse_police verbal_abuse_police		 /*
							  */ physical_abuse_military verbal_abuse_military
				
				loc j = 1
				foreach x of varlist $allabuses {
					
					if (regexm("`x'", "police") == 1) loc ti "Police"
					if (regexm("`x'", "milita") == 1) loc ti "Military"
			
					quietly reg `x' i.treatment ///
								i.barrio_code $geovars $indiv_control  [pweight=iweight], vce(cluster manzana_code)
						
						estimates store abua`j',  title("`ti'")
						estadd matrix table = r(table)
						
							loc ++j		
				}

		restore
			
			
//-----//		
//-- Coefplot of all regressions
//-----//

//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(8,3,.)
matrix S = J(8,3,.)

foreach n of numlist 1/8{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore abud1
	} 
	
	if `n' == 2 {
		estimates restore abud3
	} 
	
	if `n' == 3 {
		estimates restore abud2
	}
	
	if `n' == 4 {
		estimates restore abud4
	}
	
	if `n' == 5 {
		estimates restore abua1
	} 
	
	if `n' == 6 {
		estimates restore abua3
	} 
	
	if `n' == 7 {
		estimates restore abua2
	}
	
	if `n' == 8 {
		estimates restore abua4
	}
	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "poli_phys_dur" "mili_phys_dur" "poli_verb_dur" "mili_verb_dur" "poli_phys_aft" "mili_phys_aft" "poli_verb_aft" "mili_verb_aft"  
	matrix colname T = coef ll ul
	matrix rowname T = "poli_phys_dur" "mili_phys_dur" "poli_verb_dur" "mili_verb_dur" "poli_phys_aft" "mili_phys_aft" "poli_verb_aft" "mili_verb_aft"  

	
//-------//
//-- Run regressions for randomization inference (monitoring data - during intervention)
//-------//

* Define globals for regressions
	global indiv_control age gender 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 	 


* Run regressions

foreach var in physical_abuse_police verbal_abuse_police 		 /*
				*/ physical_abuse_military verbal_abuse_military {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
	 
	 
	//- Merge with outcomes
			merge 1:m manzana_code using "$data/raw/survey_monitoring.dta"
				drop _merge
	
	//-- Sample during the intervention (before november 19, 2019)
	
			keep if date < mdy(11,19,2019)
	
					
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${geovars} ${indiv_control} [pweight=iweight], ///
				vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA23_during_`var'.dta", replace
	
}

	
						
//-------//
//-- Run regressions for randomization inference (monitoring data - after intervention)
//-------//

* Define globals for regressions
	global indiv_control age gender 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 	 


* Run regressions

foreach var in physical_abuse_police verbal_abuse_police 		 /*
				*/ physical_abuse_military verbal_abuse_military {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
			
	 
	//- Merge with outcomes
			merge 1:m manzana_code using "$data/raw/survey_monitoring.dta"
				drop _merge
	
	//-- Sample after the intervention (after november 19, 2019)
	
			keep if date >= mdy(11,19,2019)
	
					
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${geovars} ${indiv_control} [pweight=iweight], ///
				vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA23_after_`var'.dta", replace
	
}



//------//					 
//-- Store p-values from randomization inference 
//------//

loc n = 0	
foreach var in physical_abuse_police physical_abuse_military /*
				*/ verbal_abuse_police verbal_abuse_military /*
				*/ physical_abuse_police physical_abuse_military /*
				*/ verbal_abuse_police verbal_abuse_military  {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	loc ++n
		if (`n' <= 4) {
			use "$data/rand_inf/RI_figureA23_during_`var'.dta", clear
			} 
		
		else {
			use "$data/rand_inf/RI_figureA23_after_`var'.dta", clear
			}
			
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore abud1
	} 
	
	if `n' == 2 {
		estimates restore abud3
	} 
	
	if `n' == 3 {
		estimates restore abud2
	}
	
	if `n' == 4 {
		estimates restore abud4
	}
	
	if `n' == 5 {
		estimates restore abua1
	} 
	
	if `n' == 6 {
		estimates restore abua3
	} 
	
	if `n' == 7 {
		estimates restore abua2
	}
	
	if `n' == 8 {
		estimates restore abua4
	}

	
	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f") 
	di `p_spill`n''
} 
		

//-- Export table 

	loc var1 `"\textbf{\shortstack{Police}}"'
	loc var2 `"\textbf{\shortstack{Military}}"'
	
	esttab abud1 abud3 abud2 abud4 abua1 abua3 abua2 abua4 ///
		using "$tables/table_figure_A23.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${indiv_control}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${c_geo}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var1'"' `"`var2'"' `"`var1'"' `"`var2'"' `"`var1'"' `"`var2'"' )  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{4}{c}{\textbf{During intervention}} & \multicolumn{4}{c}{\textbf{After intervention}} \\ 		///
			\cline{2-5} \cline{6-9} 					 		///
			& \multicolumn{2}{c}{\textbf{Physical}} & \multicolumn{2}{c}{\textbf{Verbal}} & \multicolumn{2}{c}{\textbf{Physical}} & \multicolumn{2}{c}{\textbf{Verbal}} \\ 		///
			\cline{2-3} \cline{4-5} \cline{6-7} \cline{8-9}  							///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )	
		
		
//-- Coefplot of all regressions 

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.25) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.25) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("poli_phys_dur" "mili_verb_dur" = `""{bf:During intervention}" "{bf:(before November 19)}""' ///
					  "poli_phys_aft" "mili_verb_aft" = `""{bf:After intervention}" "{bf:(after November 19)}""') ///
			 headings("poli_phys_dur" = "{bf:Physical abuse}" ///
					  "poli_verb_dur" = "{bf:Verbal abuse}" ///
					  "poli_phys_aft" = "{bf:Physical abuse}" ///
					  "poli_verb_aft" = "{bf:Verbal abuse}") ///
			 coeflabels("poli_phys_dur" = "Police" "mili_phys_dur" = "Military" ///
						 "poli_verb_dur" = "Police" "mili_verb_dur" = "Military" ///
						 "poli_phys_aft" = "Police" "mili_phys_aft" = "Military" ///
						 "poli_verb_aft" = "Police" "mili_verb_aft" = "Military") ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
	
graph export "$figures/figure_A23.pdf", replace	
						 
					 
					 
				//-------------------------------------------------//
				//-- 	Figure A24. Treatment effects on abuses  --//
				// 	  disaggregated by date of survey with buffer  //
				//-------------------------------------------------//
	
//-- Load data on survey data 
	
use "$data/raw/survey_monitoring.dta", clear
global indiv_control age gender 
global geovars number_buildings_sampling area bat_min cai_min ptr_min 	 
			
	
//-----//
//-- Regressions with abuses during the intervention with buffer (before December 2)
//-----//

//-- Run regressions 

		preserve
			keep if date <= mdy(12,02,2019)
			
			
				global allabuses physical_abuse_police verbal_abuse_police		 /*
							  */ physical_abuse_military verbal_abuse_military
				
				loc j = 1
				foreach x of varlist $allabuses {
					
					if (regexm("`x'", "police") == 1) loc ti "Police"
					if (regexm("`x'", "milita") == 1) loc ti "Military"
			
					quietly reg `x' i.treatment ///
								i.barrio_code $geovars $indiv_control  [pweight=iweight], vce(cluster manzana_code)
						
						estimates store abu2wd`j',  title("`ti'")
						estadd matrix table = r(table)
						
							loc ++j		
				}
			
		restore
		
	
//-----//
//-- Regressions with abuses after the intervention with buffer (after December 2)
//-----//

//-- Run regressions
	
		preserve
			keep if date > mdy(12,02,2019)
			
				global allabuses physical_abuse_police verbal_abuse_police 		 /*
							  */ physical_abuse_military verbal_abuse_military
				
				loc j = 1
				foreach x of varlist $allabuses {
					
					if (regexm("`x'", "police") == 1) loc ti "Police"
					if (regexm("`x'", "milita") == 1) loc ti "Military"
			
					quietly reg `x' i.treatment ///
								i.barrio_code $geovars $indiv_control  [pweight=iweight], vce(cluster manzana_code)
						
						estimates store abu4wd`j',  title("`ti'")
						estadd matrix table = r(table)
						
							loc ++j		
				}

		restore
	
	
			
//-----//		
//-- Coefplot of all regressions
//-----//

//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(8,3,.)
matrix S = J(8,3,.)

foreach n of numlist 1/8{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore abu2wd1
	} 
	
	if `n' == 2 {
		estimates restore abu2wd3
	} 
	
	if `n' == 3 {
		estimates restore abu2wd2
	}
	
	if `n' == 4 {
		estimates restore abu2wd4
	}
	
	if `n' == 5 {
		estimates restore abu4wd1
	} 
	
	if `n' == 6 {
		estimates restore abu4wd3
	} 
	
	if `n' == 7 {
		estimates restore abu4wd2
	}
	
	if `n' == 8 {
		estimates restore abu4wd4
	}
	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "poli_phys_dur" "mili_phys_dur" "poli_verb_dur" "mili_verb_dur" "poli_phys_aft" "mili_phys_aft" "poli_verb_aft" "mili_verb_aft"  
	matrix colname T = coef ll ul
	matrix rowname T = "poli_phys_dur" "mili_phys_dur" "poli_verb_dur" "mili_verb_dur" "poli_phys_aft" "mili_phys_aft" "poli_verb_aft" "mili_verb_aft"  


//-------//
//-- Run regressions for randomization inference (monitoring data - during intervention with buffer)
//-------//

* Define globals for regressions
	global indiv_control age gender 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 	 


* Run regressions

foreach var in physical_abuse_police verbal_abuse_police 		 /*
				*/ physical_abuse_military verbal_abuse_military {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
			
	 
	//- Merge with outcomes
			merge 1:m manzana_code using "$data/raw/survey_monitoring.dta"
				drop _merge
		
	
	//-- Sample during the intervention (before december 2, 2019)
	
			keep if date <= mdy(12,02,2019)
	
					
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${geovars} ${indiv_control} [pweight=iweight], ///
				vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA24_during_`var'.dta", replace
	
}

	
						
//-------//
//-- Run regressions for randomization inference (monitoring data - after intervention with buffer)
//-------//

* Define globals for regressions
	global indiv_control age gender 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 	 


* Run regressions

foreach var in physical_abuse_police verbal_abuse_police 		 /*
				*/ physical_abuse_military verbal_abuse_military {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
			
	 
	//- Merge with outcomes
			merge 1:m manzana_code using "$data/raw/survey_monitoring.dta"
				drop _merge
			
	
	//-- Sample after the intervention (after december 2, 2019)
	
			keep if date > mdy(12,02,2019)
	
					
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${geovars} ${indiv_control} [pweight=iweight], ///
				vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA24_after_`var'.dta", replace
	
}



	
//------//					 
//-- Store p-values from randomization inference 
//------//

loc n = 0	
foreach var in physical_abuse_police physical_abuse_military /*
				*/ verbal_abuse_police verbal_abuse_military /*
				*/ physical_abuse_police physical_abuse_military /*
				*/ verbal_abuse_police verbal_abuse_military  {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	loc ++n
		if (`n' <= 4) {
			use "$data/rand_inf/RI_figureA24_during_`var'.dta", clear
			} 
		
		else {
			use "$data/rand_inf/RI_figureA24_after_`var'.dta", clear
			}
			
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore abu2wd1
	} 
	
	if `n' == 2 {
		estimates restore abu2wd3
	} 
	
	if `n' == 3 {
		estimates restore abu2wd2
	}
	
	if `n' == 4 {
		estimates restore abu2wd4
	}
	
	if `n' == 5 {
		estimates restore abu4wd1
	} 
	
	if `n' == 6 {
		estimates restore abu4wd3
	} 
	
	if `n' == 7 {
		estimates restore abu4wd2
	}
	
	if `n' == 8 {
		estimates restore abu4wd4
	}
	
	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f") 
	di `p_spill`n''
} 
		

//-- Export table 

	loc var1 `"\textbf{\shortstack{Police}}"'
	loc var2 `"\textbf{\shortstack{Military}}"'
	
	esttab abu2wd1 abu2wd3 abu2wd2 abu2wd4 abu4wd1 abu4wd3 abu4wd2 abu4wd4 ///
		using "$tables/table_figure_A24.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${indiv_control}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${c_geo}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var1'"' `"`var2'"' `"`var1'"' `"`var2'"' `"`var1'"' `"`var2'"' )  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{4}{c}{\textbf{During intervention}} & \multicolumn{4}{c}{\textbf{After intervention}} \\ 		///
			\cline{2-5} \cline{6-9} 					 		///
			& \multicolumn{2}{c}{\textbf{Physical}} & \multicolumn{2}{c}{\textbf{Verbal}} & \multicolumn{2}{c}{\textbf{Physical}} & \multicolumn{2}{c}{\textbf{Verbal}} \\ 		///
			\cline{2-3} \cline{4-5} \cline{6-7} \cline{8-9}  							///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )	
				
		
//-- Coefplot of all regressions 

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.25) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.25) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("poli_phys_dur" "mili_verb_dur" = `""{bf:During intervention with}" "{bf:buffer (before December 2)}""' ///
					  "poli_phys_aft" "mili_verb_aft" = `""{bf:After intervention with}" "{bf:buffer (after December 2)}""') ///
			 headings("poli_phys_dur" = "{bf:Physical abuse}" ///
					  "poli_verb_dur" = "{bf:Verbal abuse}" ///
					  "poli_phys_aft" = "{bf:Physical abuse}" ///
					  "poli_verb_aft" = "{bf:Verbal abuse}") ///
			 coeflabels("poli_phys_dur" = "Police" "mili_phys_dur" = "Military" ///
						 "poli_verb_dur" = "Police" "mili_verb_dur" = "Military" ///
						 "poli_phys_aft" = "Police" "mili_phys_aft" = "Military" ///
						 "poli_verb_aft" = "Police" "mili_verb_aft" = "Military") ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
	
graph export "$figures/figure_A24.pdf", replace	
						 
						 
						 
				//-----------------------------------------//
				//--  Figure A.25 Treatment effects on   --//
				//	   exposure to police and military     //
				//-----------------------------------------//

//-----//
//-- Regressions with monitoring data (during intervention)
//-----//
						
//-- Load monitoring data  

	use "$data/raw/survey_monitoring.dta", clear
	global indiv_control age gender 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 	 

		
				
//- Seen police/military on the block 
		
	qui reg d_seen_police i.treatment ///
				i.barrio_code $geovars $indiv_control [pweight=iweight], vce(cluster manzana_code)
				
		estimates store poli_seen,  title("Police")
		estadd matrix table = r(table)
	
	qui reg d_seen_military i.treatment ///
				i.barrio_code $geovars $indiv_control [pweight=iweight], vce(cluster manzana_code)
		
		estimates store mili_seen,  title("Military")
		estadd matrix table = r(table)
	
		
//-----//
//-- Regressions with survey data (after intervention)
//-----//

//-- Load data on survey data 

	use "$data/raw/survey_endline.dta", clear

	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 


	
 //- Seen police/military on the block 
		
	qui reg d_seen_p i.treatment ///
				i.barrio_code ${demovars} ${geovars} [pweight=iweight], vce(cluster manzana_code) baselevels
		
		estimates store compl_poli,  title("Police")
		estadd matrix table = r(table)
	
	
	qui reg d_seen_m i.treatment ///
				i.barrio_code ${demovars} ${geovars} [pweight=iweight], vce(cluster manzana_code) baselevels
		
		estimates store compl_mili,  title("Military")
		estadd matrix table = r(table)
	 							
							
//-----//		
//-- Coefplot of all regressions
//-----//
 
//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(4,3,.)
matrix S = J(4,3,.)

foreach n of numlist 1/4{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore poli_seen
	} 
	
	if `n' == 2 {
		estimates restore mili_seen
	} 
	
	if `n' == 3 {
		estimates restore compl_poli
	}
	
	if `n' == 4 {
		estimates restore compl_mili
	}
	

	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "poli_dur" "mili_dur" "poli_aft" "mili_aft"   
	matrix colname T = coef ll ul
	matrix rowname T = "poli_dur" "mili_dur" "poli_aft" "mili_aft"     

							
//-------//
//-- Run regressions for randomization inference (monitoring data)
//-------//

* Define globals for regressions
	global indiv_control age gender 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 	


* Run regressions

foreach var in d_seen_police d_seen_military {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code 
			
	 
	//- Merge with outcomes
			merge 1:m manzana_code using "$data/raw/survey_monitoring.dta"
				drop _merge
			
					
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${geovars} ${indiv_control} [pweight=iweight], ///
				vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA25_`var'.dta", replace
	
}

							
						
//-------//
//-- Run regressions for randomization inference (endline data)
//-------//

* Define globals for regressions
	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 


* Run regressions

foreach var in d_seen_p d_seen_m {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
			
	 
	//- Merge with outcomes, weights and controls
		
			* Merge with endline database
			merge 1:m manzana_code using "$data/raw/survey_endline.dta"
				keep if _merge == 3
				drop _merge
	
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA25_`var'.dta", replace
	
}

						
//------//					 
//-- Store p-values from randomization inference 
//------//

loc n = 0	
foreach var in d_seen_police d_seen_military d_seen_p d_seen_m  {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	loc ++n
	use "$data/rand_inf/RI_figureA25_`var'.dta", clear	
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore poli_seen
	} 
	
	if `n' == 2 {
		estimates restore mili_seen
	} 
	
	if `n' == 3 {
		estimates restore compl_poli
	}
	
	if `n' == 4 {
		estimates restore compl_mili
	}
	
	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f") 
	di `p_spill`n''
} 
	
	
//-- Export table 

	loc var1 `"\textbf{\shortstack{Police seen on block}}"'
	loc var2 `"\textbf{\shortstack{Military seen on block}}"'
	
	esttab poli_seen mili_seen compl_poli compl_mili ///
		using "$tables/table_figure_A25.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${indiv_control}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${c_geo}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var1'"' `"`var2'"')  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{2}{c}{\textbf{During intervention}} & \multicolumn{2}{c}{\textbf{After intervention}} \\ 		///
			\cline{2-3} \cline{4-5} 					 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )
	   
	   
//-- Coefplot of all regressions 

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.08) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.08) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("poli_dur" "mili_dur" = `""{bf:During intervention}" "{bf:(monitoring survey)}""'  ///
					  "poli_aft" "mili_aft" = `""{bf:After intervention}" "{bf:(endline survey)}""') ///
			 coeflabels("poli_dur" = `""Police seen" "on block""' "mili_dur" = `""Military seen" "on block""' ///
						 "poli_aft" = `""Police seen" "on block""' "mili_aft" = `""Military seen" "on block""')  ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
	
graph export "$figures/figure_A25.pdf", replace	
						 
						 						
												
				//----------------------------------------------//
				//-- Figure A.26 Treatment effects on police  --//
				//	and military activities during intervention	//
				//----------------------------------------------//
				
//-----//
//-- Regressions with monitoring data (during intervention)
//-----//
						
//-- Load monitoring data  

	use "$data/raw/survey_monitoring.dta", clear
	global indiv_control age gender 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 
	
	
//- Seen police/military making arrests 
	
	qui reg d_arresting_police i.treatment 	///
			i.barrio_code $geovars $indiv_control [pweight = iweight], vce(cluster manzana_code)
		
		estimates store poli_arrest, title("Police")
		estadd matrix table = r(table)
	
	
	qui reg d_arresting_military i.treatment 	///
			i.barrio_code $geovars $indiv_control [pweight = iweight], vce(cluster manzana_code)
		
		estimates store mili_arrest, title("Military")
		estadd matrix table = r(table)
			   
			
//- Seen police/military checking ids 
			
	qui reg d_idchecking_police i.treatment 	///
			i.barrio_code $geovars $indiv_control [pweight = iweight], vce(cluster manzana_code)
		
		estimates store poli_ch, title("Police")
		estadd matrix table = r(table)
		
	
	qui reg d_idchecking_military i.treatment 	///
			i.barrio_code $geovars $indiv_control [pweight = iweight], vce(cluster manzana_code)
		
		estimates store mili_ch, title("Military")
		estadd matrix table = r(table)
	
	
//- Seen police/miltiary talking with citizens 
	
	qui reg d_interacting_police  i.treatment 	///
			i.barrio_code $geovars $indiv_control [pweight = iweight], vce(cluster manzana_code)
		
		estimates store poli_int, title("Police")
		estadd matrix table = r(table)
		
	
	qui reg d_interacting_military  i.treatment 	///
			i.barrio_code $geovars $indiv_control [pweight = iweight], vce(cluster manzana_code)
		
		estimates store mili_int, title("Military")
		estadd matrix table = r(table)

		
//-----//		
//-- Coefplot of all regressions
//-----//

//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(6,3,.)
matrix S = J(6,3,.)

foreach n of numlist 1/6{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore poli_arrest
	} 
	
	if `n' == 2 {
		estimates restore mili_arrest
	} 
	
	if `n' == 3 {
		estimates restore poli_ch
	}
	
	if `n' == 4 {
		estimates restore mili_ch
	}
	
	if `n' == 5 {
		estimates restore poli_int
	} 
	
	if `n' == 6 {
		estimates restore mili_int
	} 
	
	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "poli_arrest" "mili_arrest" "poli_ch" "mili_ch" "poli_int" "mili_int" 
	matrix colname T = coef ll ul
	matrix rowname T = "poli_arrest" "mili_arrest" "poli_ch" "mili_ch" "poli_int" "mili_int" 

			
//-------//
//-- Run regressions for randomization inference (monitoring data)
//-------//

* Define globals for regressions
	global indiv_control age gender 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 


* Run regressions

foreach var in d_arresting_police d_arresting_military d_idchecking_police d_idchecking_military /*
			*/ d_interacting_police d_interacting_military {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code 
			
	 
	//- Merge with outcomes
			merge 1:m manzana_code using "$data/raw/survey_monitoring.dta"
				keep if _merge == 3
				drop _merge
			
					
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${geovars} ${indiv_control} [pweight=iweight], ///
				vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA26_during_`var'.dta", replace
	
}



//------//					 
//-- Store p-values from randomization inference 
//------//

loc n = 0	
foreach var in d_arresting_police d_arresting_military d_idchecking_police d_idchecking_military /*
			*/ d_interacting_police d_interacting_military {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	loc ++n
	use "$data/rand_inf/RI_figureA26_during_`var'.dta", clear
			
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore poli_arrest
	} 
	
	if `n' == 2 {
		estimates restore mili_arrest
	} 
	
	if `n' == 3 {
		estimates restore poli_ch
	}
	
	if `n' == 4 {
		estimates restore mili_ch
	}
	
	if `n' == 5 {
		estimates restore poli_int
	} 
	
	if `n' == 6 {
		estimates restore mili_int
	} 
	
	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f") 
	di `p_spill`n''
} 
	
	
//-- Export table 

	loc var1 `"\textbf{\shortstack{Police}}"'
	loc var2 `"\textbf{\shortstack{Military}}"'
	
	esttab poli_arrest mili_arrest poli_ch mili_ch poli_int mili_int ///
		using "$tables/table_figure_A26.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${indiv_control}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${c_geo}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var1'"' `"`var2'"'  `"`var1'"' `"`var2'"')  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{6}{c}{\textbf{During intervention}} \\ 		///
			\cline{2-7} 								 		///
			& \multicolumn{2}{c}{\textbf{Seen making arrests}} & \multicolumn{2}{c}{\textbf{Seen checking ID}} & \multicolumn{2}{c}{\textbf{Seen talking with citizens}}  \\ 		///
			\cline{2-3} \cline{4-5} \cline{6-7} 				///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )
	   
	  
	  
//-- Coefplot of all regressions 

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.17) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.17) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xscale(range(-0.05 0.14)) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("poli_arrest" "mili_int" = `""{bf:During intervention}" "{bf:(monitoring survey)}""') ///
			 headings("poli_arrest" = "{bf:Seen making arrests}" ///
					  "poli_ch" = "{bf:Seen checking ID}" ///
					  "poli_int" = "{bf:Seen talking with citizens}") ///
			 coeflabels("poli_arrest" = "Police" "mili_arrest" = "Military" ///
						 "poli_ch" = "Police" "mili_ch" = "Military" ///
						 "poli_int" = "Police" "mili_int" = "Military" )  ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
	
graph export "$figures/figure_A26.pdf", replace	
						 		
			

				//----------------------------------------------//
				//--  Figure A.27 Treatment effects on police --//
				//	and military activities after intervention	//
				//----------------------------------------------//

//-----//
//-- Regressions with survey data (after intervention)
//-----//

//-- Load data on survey data 

	use "$data/raw/survey_endline.dta", clear

	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 


//- Seen police/military making arrests 
		
	qui reg d_arresting_p i.treatment ///
				i.barrio_code ${demovars} ${geovars} [pweight=iweight], vce(cluster manzana_code) baselevels
		
		estimates store poli_arrest_aft,  title("Police")
		estadd matrix table = r(table)
	
	
	qui reg d_arresting_m i.treatment ///
				i.barrio_code ${demovars} ${geovars} [pweight=iweight], vce(cluster manzana_code) baselevels
		
		estimates store mili_arrest_aft,  title("Military")
		estadd matrix table = r(table)
	
	
//- Seen police/military checking IDs 
		
	qui reg d_idchecking_p i.treatment ///
				i.barrio_code ${demovars} ${geovars} [pweight=iweight], vce(cluster manzana_code) baselevels
		
		estimates store poli_ch_aft,  title("Police")
		estadd matrix table = r(table)
	
	
	qui reg d_idchecking_m i.treatment ///
				i.barrio_code ${demovars} ${geovars} [pweight=iweight], vce(cluster manzana_code) baselevels
		
		estimates store mili_ch_aft,  title("Military")
		estadd matrix table = r(table)
	
	
//- Seen police/military talking with citizens
		
	qui reg d_interacting_p i.treatment ///
				i.barrio_code ${demovars} ${geovars} [pweight=iweight], vce(cluster manzana_code) baselevels
		
		estimates store poli_int_aft,  title("Police")
		estadd matrix table = r(table)
	
	
	qui reg d_interacting_m i.treatment ///
				i.barrio_code ${demovars} ${geovars} [pweight=iweight], vce(cluster manzana_code) baselevels
		
		estimates store mili_int_aft,  title("Military")
		estadd matrix table = r(table)
			
	
//-----//		
//-- Coefplot of all regressions
//-----//

//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(6,3,.)
matrix S = J(6,3,.)

foreach n of numlist 1/6{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore poli_arrest_aft
	}
	
	if `n' == 2 {
		estimates restore mili_arrest_aft
	}
	
	if `n' == 3 {
		estimates restore poli_ch_aft
	}
	
	if `n' == 4 {
		estimates restore mili_ch_aft
	} 
	
	if `n' == 5 {
		estimates restore poli_int_aft
	} 
	
	if `n' == 6 {
		estimates restore mili_int_aft
	}
	
	
	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "poli_arrest_aft" "mili_arrest_aft" "poli_ch_aft" "mili_ch_aft" "poli_int_aft" "mili_int_aft"
	matrix colname T = coef ll ul
	matrix rowname T = "poli_arrest_aft" "mili_arrest_aft" "poli_ch_aft" "mili_ch_aft" "poli_int_aft" "mili_int_aft"

			
//-------//
//-- Run regressions for randomization inference (endline data)
//-------//

* Define globals for regressions
	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 


* Run regressions

foreach var in d_arresting_p d_arresting_m d_idchecking_p d_idchecking_m /*
			*/ d_interacting_p d_interacting_m {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
			
	 
	//- Merge with outcomes, weights and controls
		
			* Merge with endline database
			merge 1:m manzana_code using "$data/raw/survey_endline.dta"
				keep if _merge == 3
				drop _merge
	
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA27_after_`var'.dta", replace
	
}

						
//------//					 
//-- Store p-values from randomization inference 
//------//

loc n = 0	
foreach var in d_arresting_p d_arresting_m d_idchecking_p d_idchecking_m /*
			*/ d_interacting_p d_interacting_m {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	loc ++n
	use "$data/rand_inf/RI_figureA27_after_`var'.dta", clear
			
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore poli_arrest_aft
	}
	
	if `n' == 2 {
		estimates restore mili_arrest_aft
	}
	
	if `n' == 3 {
		estimates restore poli_ch_aft
	}
	
	if `n' == 4 {
		estimates restore mili_ch_aft
	} 
	
	if `n' == 5 {
		estimates restore poli_int_aft
	} 
	
	if `n' == 6 {
		estimates restore mili_int_aft
	}
	
	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f") 
	di `p_spill`n''
} 
		

//-- Export table 

	loc var1 `"\textbf{\shortstack{Police}}"'
	loc var2 `"\textbf{\shortstack{Military}}"'
	
	esttab poli_arrest_aft mili_arrest_aft poli_ch_aft mili_ch_aft poli_int_aft mili_int_aft ///
		using "$tables/table_figure_A27.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${indiv_control}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${c_geo}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var1'"' `"`var2'"'  `"`var1'"' `"`var2'"')  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{6}{c}{\textbf{After intervention}} \\ 		///
			\cline{2-7} 								 		///
			& \multicolumn{2}{c}{\textbf{Seen making arrests}} & \multicolumn{2}{c}{\textbf{Seen checking ID}} & \multicolumn{2}{c}{\textbf{Seen talking with citizens}}  \\ 		///
			\cline{2-3} \cline{4-5} \cline{6-7} 				///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )
	  		
		
//-- Coefplot of all regressions 

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.17) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.17) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xscale(range(-0.05 0.14)) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 groups("poli_arrest_aft" "mili_int_aft" = `""{bf:After intervention}" "{bf:(endline survey)}""') ///
			 headings("poli_arrest_aft" = "{bf:Seen making arrests}" ///
					  "poli_ch_aft" = "{bf:Seen checking ID}" ///
					  "poli_int_aft" = "{bf:Seen talking with citizens}" ) ///
			 coeflabels("poli_arrest_aft" = "Police" "mili_arrest_aft" = "Military" ///
						 "poli_ch_aft" = "Police" "mili_ch_aft" = "Military" ///
						 "poli_int_aft" = "Police" "mili_int_aft" = "Military" )  ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
	
graph export "$figures/figure_A27.pdf", replace	
			

				
					//--------------------------------------------//
					//-- 	Figure A.28. Treatment effects on   --//
					// 	  cooperation with police and military    //
					//--------------------------------------------//

//-----//
//-- Regressions with survey data (after intervention)
//-----//

//-- Load data on survey data 

	use "$data/raw/survey_endline.dta", clear

	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 


//-- Run regressions

	qui reg i_coopallindex_std i.treatment ///
				i.barrio_code ${demovars} ${geovars} [pweight=iweight], vce(cluster manzana_code) baselevels
		
		estimates store coop_e,  title("All")
		estadd matrix table = r(table)
		
		
	qui reg i_cooppoliceindex_std i.treatment ///
				i.barrio_code ${demovars} ${geovars} [pweight=iweight], vce(cluster manzana_code) baselevels
		
		estimates store coop_pe,  title("Police")
		estadd matrix table = r(table)
		
	
	qui reg i_coopmilitaryindex_std i.treatment ///
				i.barrio_code ${demovars} ${geovars} [pweight=iweight], vce(cluster manzana_code) baselevels
		
		estimates store coop_me,  title("Military")	
		estadd matrix table = r(table)
	

//-----//		
//-- Coefplot of all regressions
//-----//

//-- Obtain coefficients and confidence intervals of each regression 

matrix T = J(3,3,.)
matrix S = J(3,3,.)

foreach n of numlist 1/3{
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore coop_e
	} 
	
	if `n' == 2 {
		estimates restore coop_pe
	} 
	
	if `n' == 3 {
		estimates restore coop_me
	}
	
	
	* Store estimates in matrix 
	
	matrix t_table = e(table)
			matrix T[`n',1] = t_table[1,"1.treatment"]
			matrix T[`n',2] = t_table[5,"1.treatment"]
			matrix T[`n',3] = t_table[6,"1.treatment"]
			
			matrix S[`n',1] = t_table[1,"2.treatment"]
			matrix S[`n',2] = t_table[5,"2.treatment"]
			matrix S[`n',3] = t_table[6,"2.treatment"]
}


	* Name matrix elements 
	
	matrix colname S = coef ll ul
	matrix rowname S = "all" "police" "military"
	matrix colname T = coef ll ul
	matrix rowname T = "all" "police" "military"
	
							
//-------//
//-- Run regressions for randomization inference (endline data)
//-------//

* Define globals for regressions
	global demovars age gender educ 
	global geovars number_buildings_sampling area bat_min cai_min ptr_min 


* Run regressions

foreach var in i_coopallindex_std i_cooppoliceindex_std i_coopmilitaryindex_std {
	
	di "*** REGRESSIONS WITH VARIABLE `var' ***"
	
	* Create empty matrix
	matrix treat_effects = J(10000, 2, .)
	matrix colnames treat_effects = treat_ef spillover_ef
	
	
	forvalues n = 1/10000{
	
			use "$data/rand_inf/block_simulate_randomizations_p1.dta", clear
				
			di "This is simulation number `n'"
			
		
	//- Obtain treatment status of each block 
			keep manzana treatment_ri_`n'
			rename treatment_ri_`n' treatment_ri
			rename manzana manzana_code
			
	 
	//- Merge with outcomes
			merge 1:m manzana_code using "$data/raw/survey_endline.dta"
				drop _merge
		
			
	//- Run regressions 
			qui: reg `var' i.treatment_ri i.barrio_code ${demovars} ${geovars} [pweight=iweight], ///
			vce(cluster manzana_code) baselevels
	
			
			* Store coefficients in matrix
			matrix treat_effects[`n', 1] = _b[1.treatment_ri]
			matrix treat_effects[`n', 2] = _b[2.treatment_ri]
}

* Turn ATEs into variables
	
	clear	
	
	svmat treat_effects, names(col)
	
	save "$data/rand_inf/RI_figureA28_`var'.dta", replace
	
}

						

//------//					 
//-- Store p-values from randomization inference 
//------//

loc n = 0	
foreach var in i_coopallindex_std i_cooppoliceindex_std i_coopmilitaryindex_std {
	
	* Load data of simulations
	
	di "*** THIS IS COLUMN NUMBER `var' ***"
	
	loc ++n
	use "$data/rand_inf/RI_figureA28_`var'.dta", clear
			
	
	* Restore estimates of each regression
	
	if `n' == 1 {
		estimates restore coop_e
	} 
	
	if `n' == 2 {
		estimates restore coop_pe
	} 
	
	if `n' == 3 {
		estimates restore coop_me
	}
	
	* P-value of treatment efffect

	count if abs(treat_ef) > abs(_b[1.treatment])
	estadd scalar p_val_treat =  `r(N)'/10000, replace
	
	local p_treat`n' = string(`e(p_val_treat)', "%9.3f")
	di `p_treat`n'' 
	
	
	* P-value of spillover effect

	count if abs(spillover_ef) > abs(_b[2.treatment])
	estadd scalar p_val_spillover =  `r(N)'/10000, replace
	
	local p_spill`n' = string(`e(p_val_spillover)', "%9.3f") 
	di `p_spill`n''
} 
		
		
//-- Export table 

	loc var1 `"\textbf{\shortstack{All}}"'
	loc var2 `"\textbf{\shortstack{Police}}"'
	loc var3 `"\textbf{\shortstack{Military}}"'
	
	esttab coop_e coop_pe coop_me ///
		using "$tables/table_figure_A28.tex", ///
		legend booktabs f replace /// ///
		order(1.treatment 2.treatment) ///
		keep(1.treatment 2.treatment) /// 
		indicate("Individual controls = ${indiv_control}" "Neighborhood FE = *.barrio_code*" "Block-level controls = ${c_geo}" , labels("Yes" "No")) ///
		stats(N r2 p_val_treat p_val_spillover, fmt(0 2 3 3 0 0)  labels("Observations" "\$R^2$" "RI p-value (treatment)" "RI p-value (spillover)")) ///
		nostar	///
		cells("b(fmt(3))" "ci_l(par) & ci_u(par)" "p(fmt(3) par)") incelldelimiter(", ") ///
		label nonotes nonum nodepvar nogaps mtitle(`"`var1'"' `"`var2'"' `"`var3'"')  ///
		collabels(none) ///
		nolines														///
		prehead(  												///
		   \def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}			///
		   \begin{tabular}{@{\extracolsep{4pt}}l*{@M}{c}@{}} 	/// 
		   \hline \hline     									///
		   \noalign{\smallskip}									///
		   & \multicolumn{3}{c}{\textbf{Cooperation with authorities}} \\ 		///
			\cline{2-4} 								 		///
	   )  														///
	   posthead(\addlinespace 									///
			\hline 												///
			\addlinespace										///
		) 														///
	   postfoot(  												///
		  \noalign{\smallskip} \hline \hline  					///
		  \end{tabular}											///
		  \medskip    											///
	   )
				
				
//-- Coefplot of all regressions 

coefplot (matrix(T[.,1]), ci((T[.,2] T[.,3])) label(Treatment) offset(0.05) pstyle(p3) msymbol(O) mcolor(black) ciopts(lpatt(solid) lcol(black)))   ///
			 (matrix(S[.,1]), ci((S[.,2] S[.,3])) label(Spillover) offset(-0.05) pstyle(p4) msymbol(S) mcolor(gs5) ciopts(lpatt(dash) lcol(gs5))), ///
			 graphregion(fcolor(white)) ///
			 ciopts(recast(rcap)) ///
			 xline(0, lp(dot) lc(black)) ///
			 scheme(plotplain) legend(position(3) rows(2)) ///
			 xlabel(, nogrid) xscale(range(-0.05 0.25)) xtitle("ITT effect size", size(small)) ///
			 ylabel(, nogrid) ///
			 coeflabels("all" = "All" "police" = "Police" "military" = "Military") 	///
			 groups("all" "military" = `""{bf:Cooperation with authorities}" "{bf:after intervention}" "{bf:(endline survey)}""')  ///
			mlabcolor(none) addplot(scatter @at @ul, ms(i) mlabel(@mlbl) mlabcolor(black))
	
graph export "$figures/figure_A28.pdf", replace	
			
									
						
						